/**
 * @file   Matrix.h
 * @author Wang Heyu <wang@wukong>
 * @date   Tue May  4 21:12:37 2021
 *
 * @brief  Dense Matrix.
 *
 *
 */

#ifndef __CRAZYFISH_MATRIX__
#define __CRAZYFISH_MATRIX__

#include <valarray>
#include <cstdlib>
#include <iostream>

#define TEMPLATE template <typename T>

TEMPLATE class Matrix;

TEMPLATE
class Matrix : public std::valarray<T>
{
private:
    size_t n_col;

public:
    /**
     * Constructor. Of giving the row and column sizes by @p _m and @p
     * _n , respectively.
     *
     * @param _m row size.
     * @param _n column size.
     *
     */
    Matrix(size_t _m, size_t _n);
    T &operator()(size_t _i, size_t _j);
    const T &operator()(size_t _i, size_t _j) const;
    std::valarray<T> operator*(const std::valarray<T> &_p) const;
    size_t get_n_row() const;
    size_t get_n_col() const;
    Matrix<T> transpose() const;
};

TEMPLATE
Matrix<T>::Matrix(size_t _m, size_t _n) : std::valarray<T>(_m * _n), // 若T类型不为数字，则会报错
                                          n_col(_n){
                                              //    n_col = _n;
                                          };

TEMPLATE
size_t Matrix<T>::get_n_col() const
{
    return n_col;
};

TEMPLATE
size_t Matrix<T>::get_n_row() const
{
    if (n_col == 0)
    {
        return 0;
    }
    else
    {
        return (this->size() / n_col);
    }
};

TEMPLATE
T &Matrix<T>::operator()(size_t _i, size_t _j)
{
    if ((_j >= 0) && (_j < n_col))
    {
        return (*this)[_i * n_col + _j];
    }
    else
    {
        std::cerr << "index must be positive and less the scale of Matrix." << std::endl;
        exit(-1);
    }
};

TEMPLATE
const T &Matrix<T>::operator()(size_t _i, size_t _j) const
{
    if ((_j >= 0) && (_j < n_col))
    {
        return (*this)[_i * n_col + _j];
    }
    else
    {
        std::cerr << "index must be positive and less the scale of Matrix." << std::endl;
        exit(-1);
    }
};

TEMPLATE
std::valarray<T> Matrix<T>::operator*(const std::valarray<T> &_p) const
{
    size_t n_row = get_n_row();
    size_t n_col = get_n_col();
    if (n_col != _p.size())
    {
        std::cerr << "Dimension between the matrix and vactor are not matching." << std::endl;
        exit(-1);
    }
    std::valarray<T> re(n_row);

    for (size_t i = 0; i < n_row; i++)
    {
        std::slice myslice = std::slice(i * n_col, n_col, 1);
        re[i] = ((*this)[myslice] * _p).sum();
    }
    return re;
};

TEMPLATE
Matrix<T> Matrix<T>::transpose() const
{
    size_t n_row = get_n_row();
    Matrix<T> T_Mat(n_row, n_col);
    for (size_t i = 0; i < n_row; i++)
        for (size_t j = 0; j < n_col; j++)
            T_Mat(i, j) = (*this)[j * n_col + i];
    return T_Mat;
}

#undef TEMPLATE
#else
// DO NOTHING.
#endif
